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Numerical study of boundary layer 
interaction with shocks — method 
improvement and test computation 


By N. A. Adams 


1. Motivation and objectives 

The general motivation of this work has been outlined in Adams (1994). The 
objective is the development of a high-order and high-resolution method for the 
direct numerical simulation of shock turbulent-boundaxy-layer interaction. Details 
concerning the spatial discretization of the convective terms can be found in Adams 
and Shariff (1995). The computer code based on this method as introduced in 
Adams (1994) was formulated in Cartesian coordinates and thus has been limited 
to simple rectangular domains. For more general two-dimensional geometries, as a 
compression corner, an extension to generalized coordinates is necessary. To keep 
the requirements or limitations for grid generation low, the extended formulation 
should allow for non- orthogonal grids. Still, for simplicity and cost efficiency, peri- 
odicity can be assumed in one cross-flow direction. 

For easy vectorization, the compact-ENO coupling algorithm as used in Adams 
(1994) treated whole planes normal to the derivative direction with the ENO scheme 
whenever at least one point of this plane satisfied the detection criterion. This is 
apparently too restrictive for more general geometries and more complex shock 
patterns. Here we introduce a localized compact-ENO coupling algorithm, which is 
efficient as long as the overall number of grid points treated by the ENO scheme is 
small compared to the total number of grid points. 

Validation and test computations with the fined code axe performed to assess the 
efficiency and suitability of the computer code for the problems of interest. We de- 
fine a set of parameters where a direct numerical simulation of a turbulent boundary 
layer along a compression corner with reasonably fine resolution is affordable. 


2. Accomplishments 


2.1 Generalized coordinates 

The fundamental equations solved are the conservation equations for mass, mo- 
mentum, and energy in generalized coordinates 

d_U_ 3 _Fe d H e _ d F s dGs d_Hs 

dt J dx J ^ dy J dz J dx J dy J dz J 


( 1 ) 
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where the conservative variables are 


U = 


P 

pu 

pv 

pw 

IE] 


( 2 ) 


with E = ^yp+|(u 2 + t; 2 +ti; 2 ). Considering only spanwise periodic configurations 
we limit the coordinate generalization to the (x, z)-plane. The convective fluxes are 
given by 

p(u£ x + w£z) 
pu(u£ x + w(z) + pi x 
pv(u(, x + w£ z ) 

pw{uZx + wZz) 

L (E + p)(u£ * + w£z) 

and similarly for Ge and He- The viscous fluxes are given by 

T xzt,z "I" Txz(>z 


F e = 


(3) 


F s = 




Tzz£z "b Tzz£ z 

■ -<htz - qziz + ( UTxz + VTxy + + ( UTxz + VT yz + WTzz)Zz j 

and similarly Gs and Hs- The Jacobian of the coordinate transformation is 

J — £z£z ~ £zGz ■ 

The stresses are defined as 


p 

'4 

f du du \ 

2 dv 

2 i 

fdw dw V 

Txx = ~Re 

3 ^ 

{W* + ac c v 

"3 d~A~ 

3< 

[dC x+ ac c v. 


with analogous definitions for T yy and t zz ; 

p\(dv dv \ du 

Tz * ~ Te + 5C C V dv Vy \ ’ 

and similarly for r IZ , r yz , and t zz . The heat fluxes are defined as 


q z = 


(?L e +?LA 

\dC x+ d( Qx ) ’ 


(4) 


( 5 ) 


( 6 ) 


(7) 


( 8 ) 


(k - l)M^PrRe 

q y and q z analogously. The viscosity is calculated according to Sutherland’s law. 
We also assume the thermal equation of state for perfect gases to be valid. 

Given a wall-normal temperature gradient distribution <9T/3n, a von Neumann 
condition for the temperature is imposed by setting 


v'EFkFK-k, (.+«.)£ 

‘ C2 + C? 


(9) 


whenever it appears during the computation of heat flux and stress terms (due to 
the temperature dependence of the viscosity). 
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2.2 Grid generation 

For the generation of an analytic mapping of the computational domain onto the 
physical domain we follow a simple algebraic procedure. We restrict our interest to 
channel-like geometries where lower and upper boundary can be approximated by 
simple functions. The mapping is non-conformal and thus the orthogonal partition 
of the computational domain will be mapped onto a non-orthogonal partition of the 
physical domain in general. The mapping consists of two steps: (1) the computa- 
tional domain {£, (} £ [0, 1] x [0, 1] with a uniformly spaced partitioning is mapped 
onto an intermediate space with non-uniform partitioning {$,£} € [0, 1] x [0, 1]; (2) 
the intermediate space {s,t} is mapped onto the physical space {x,z}. Using a 
linear blending function between lower and upper boundary, we define this latter 
mapping function by 

x(£, C) = (1 - + tx u (s) (10) 

*(£, C) = M 5 ) + (! - *K( 5 ) > ( 11 ) 


the indices / and u indicate that the functions are to be taken at the lower and 
upper boundary, respectively. The components of the Jacobi matrix are then given 

by 

d(x,z) 


d(U) 


r 

dz_ 

L a* 


dx 

dz , 


( 12 ) 


Later the metric coefficients will be needed, which are the components of the inverse 
Jacobi matrix, 


d( 


_ i j. ( f n , ( d{x,z) Y 

x,z) J U«,oJ[ v«K.<))J 


-1 


and the Jacobian 


J(t, C) = Det 


( %M) \ 

\d(x,z)J • 


(13) 

(14) 


For the point distributions along the parameter lines s(£) along the lower and 
upper boundary, we define 


and its derivative 


s(0 = a£ + 6 + ci sinh[ 0 3 (Ol 

ds(0 


^1 

= a + — cosh[jf 3 (0] • 

a? c 3 


The following abbreviations axe used: 


9^(0 = 


£ “ c 2 


C3 


a = 1 — ci 


sinh ( - 

) + smh ( ) 

\c 

3 / l Q /J 


(15) 

(16) 

(17) 

(18) 
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fc = ci sinh • (19) 

If we consider compression corner geometries, then ci and c 3 are parameters that 
tune the grid point condensation around the corner point x c . It coincides with the 
zero of sinh[< 73 (£)], which is the condition from which c 2 is computed by solving 

x c - x(c 2 ) = 0 (20) 

for c 2 . Knowing all parameters we define the variation of x along the lower or upper 
boundary in terms of the parameter s as 

x(s) = Ls , (21) 

where L is the maximum value x assumes on the lower or upper boundary, respec- 
tively. Having obtained x(s) we get z(s) in the following manner. 


z(x) = d 2 


x + — ln(cosh(dix — x c )) + d 3 
di 


( 22 ) 


A corner singularity in the mapping is avoided by prescribing a finite curvature r c 
at {x c ,0}. The ramp endpoint is given by {L,sin(<f>)(L — x c )}, where <f> is the ramp 
angle in physical space {x,z}. The parameter d 2 is computed from the condition 


Finally one sets 


and 


z(L) = sin(<£)(L — x c ) 

d (l+dpj 


d 3 = — -pln(cosh(dix c )) . 


In the transversal direction we introduce the parameter function t(() 

«>= if- wo [(>- if) -or 


and its derivative 


dm 

d C 


= if( L ir 1 )( c+ l CO8h( ' ,2(0) )' < 


(23) 

(24) 

(25) 

(26) 


(27) 

(28) 
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Herein following abbreviations and parameters are used: e\ and e 3 control the grid 
stretching at a point {0, z mv } similar to d\ and d 3 ; about half of the grid points are 
between {0,0} and {0,zi/ 2 }. {0,zj} is the upper-left corner point. The auxiliary 
functions h\ and h -2 are defined as 


h 1 (0 = c( + d + e 1 smh(h 2 (0) (29) 

and 

h 2 ( 0 = • (30) 

^3 

The constants c and d are given by 


and 


c = 1 — t\ 


5inh (V i ) +sinh (S)] 

d — t\ sinh ( — 



(31) 


(32) 


Given z mv , the parameter e 2 is computed from the condition that the argument of 
z(Q = Zmv coincides with / 12 (C) = 0, i.e. e 2 is obtained by solving 


z mv — 0 • 


(33) 


2.3 Local compact- ENO coupling 

The principle of the coupling between ENO-scheme and the compact finite- 
difference scheme is discussed in Adams and Shariff (1995). The actual imple- 
mentation with a reasonable capability for vector optimization is more involved. 
Let us consider the one dimensional and one component problem. Given the flux 
F on the grid {zj}, its derivative for x is approximated by 

( ^- = P N [F] = Ml 1 M R F (34) 

Assume that {x } e = {x p , . . . , x 9 }i U ■ • • U {x r , . . . , x,} ne is the union of regions of 
points where the flux derivatives are approximated by the ENO scheme. If a shock 
detection algorithm has detected a point x* to be treated by the ENO scheme, V{ 
is set true and we define a topology vector T by 

T = {v t } . (35) 


This vector has n e unity blocks with dimensions N„ E > 2N, ep +l, where N sep is the 
dimension of the padding on both sides of ENO regions (Adams and Shariff, 1995). 
Whenever we have = 1 for a certain grid point, we calculate P^[F]i = P^ NO [F], 
from the ENO scheme. 
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The effect of the compact-ENO coupling on Eq. (34) is that the ith component 
of M rF is replaced by the flux derivative at i calculated with the ENO scheme 
whenever V{ is 1. The row i of Ml hats then to become unity so that the ENO 
flux derivative P$ NO [F] is exactly returned when Eq. (34) is solved for Pn[F ]. 
We defining a correction matrix BCD, which changes the rows i of Ml to unity 
whenever v t is true by its dyadic decomposition into the matrices B, D, and C, the 
dimensions of which are given below. With this definition the fundamental equation 
for the computation of flux derivatives of the hybrid scheme can be written as 

(M l - BDC)P N [F] = M R F + T(P£ wo [F] - M R F) . (36) 

The rank of the correction matrix BCD is = the- It is evident that 

Eq. (36) returns the ENO flux derivatives exactly at points i whenever v x = 1. 

To solve Eq. (36) efficiently we make use of the identity by Frobenius and Schur 
(Zurmiihl and Falk, 1984, pg. 308,312) which allows to compute (Ml — BDC)* -1 
by using the inverse of Ml corrected by a the inverse of a rank hie matrix R. If 
rriE « N this procedure is more efficient for multi-dimensional problems by using 
the precomputed inverse of Ml than inverting the LHS-matrix of Eq. (36). 

The matrices B, D and C are defined as follows: 



m e 


B 

ii 

M 

* H 

(37) 

m e xN 

u= 1 


^ = 

^1^, — Im e 

(38) 

Trig x m e 

m.£ XrriE 


= 

mg xN 

B t (Ml - 1) • 

(39) 


Here we define e„ as the m#- component vector with its v component equal to unity 
the rest being zero. 

The solution algorithm for Eq. (36) according to (Zurmiihl and Falk, 1984) is the 
following: 

(0.) calculate the uncorrected solution vector y from 

M L y = M r F + T(Pf NO [F] - M r F) 

by direct inversion using the precomputed LU-decomposition of Ml; 

(1.) compute matrix V from 

MlV = B 

by direct inversion using the precomputed LU-decomposition of Ml ; 

(2.) generate the rank ms correction matrix R from 


R = I mE - CV ; 
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TABLE 1. Flow parameters for Moo = 6 ramp. 


quantity 

value 

comment 


57.32K 

free stream temperature 

Moo 

6 

free stream Mach number 

P*o o 

681.15Pa 

free stream pressure 

Pr 

0.7 

Prandtl number 

K 

1.4 

ratio of specific heats 

R 

287.03 

gas constant 

Pic, 

3.77 • 10~ 6 kg/m/s 

free stream viscosity 

S* 

110.4K 

Sutherland constant 

R e to 

100000 

running length Reynolds number 

Reft 

5255 

reference Reynolds number 


5.2554 • 10 -4 m 

reference length 

to 

19.03 

inflow dist. from lead, edge 

L i 

57.14 

length of first ramp segment 

l 2 

120 

length of second ramp segment 

<t> 

7.5° 

ramp deflection angle 


(3.) get the solution correction vector z from 

Rz = Cy 

(note that R is usually fully occupied so that this procedure is only efficient if 

tub « N); 

(4.) find the solution vector from 


fV[-F] = y + Z . 

For a multidimensional problem all points in index planes normal to the derivative 
direction are gathered and a vector loop is spanned over these. 

2.4 Code validation 

Similar to Adams (1994) we validate the generalized coordinate code by compar- 
ison with a steady state solution. Experimental and numerical data for comparison 
are taken from the computational and experimental results of a laminar boundary 
layer along a 7.5° compression corner at M 0 Q = 6 by Simeonides et aL (1994). 
We emphasize that for the results presented in this section time-accurate and low- 
dissipation methods have been used. The computations have thus been halted before 
a true steady state has been reached (residual about 10“ 4 ). The flow parameters 
are given in Table 1 (reference length is 6*, dimensional quantities axe marked with 
a star). 

In Fig. 2 the grid generated by the algorithm in section 2.2 is shown (each 4th 
grid line). The grid is condensed towards the wall and towards the kink of the 
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FIGURE 1. Skin friction and surface pressure, 7.5° laminar compression corner at 

Moo = 6. Symbols: , EN03TVDR3; , CUHDE4R3; o , Simeonides et 

al, exp.; o , Simeonides et a/., comput. 

ramp. As initial condition we take outside of the boundary layer the solution of the 
inviscid deflection problem, while near the wall a boundary layer from a similarity 
solution is given (ignoring the adverse pressure gradient on the inclined segment 
of the ramp). As boundary conditions we fix at the inflow the initial condition 
for all primitive variables giving the correct number of 5 conditions for the Navier- 
Stokes equations (Oliger & Sundstrom, 1978). At the outflow we prescribe perfectly 
non-reflecting boundary conditions (Thompson, 1987). At the upper boundary 
freest ream conditions for all flow variables are prescribed. 

The computation is started with N x = 151 and N z = 61. After 1000 iterations 
with a 3rd order LLF-ENO scheme, the resolution is increased to N x = 351 and 
N z = 121 and the computation is continued for 12000 time steps. Finally, we switch 
to the hybrid scheme (5th compact upwind, 4th order LLF-ENO) and continue for 
another 16000 iterations. For the shock detection parameters we use /3 X = 5 and 
f} z = 5. The agreement between the computational and experimental results of 
Simeonides et al. (1994) and the present results is satisfactory, Fig. 1. A small 
inflow transient is caused by the fact that we prescribe a boundary layer profile at 
inflow. This is to match the procedure in later DNS. In Simeonides et al. (1994) 
the plate leading edge is included in the computational domain. 

Figure 3 shows a quasi- Schlieren plot (merely the norm of the density gradient) 
when the computations were halted. Both the separation shock and the main ramp 
shock are clearly visible. 
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FIGURE 2. Grid for 7.5° ramp, each 4th grid line shown. 



FIGURE 3. Quasi- Schlieren plot (intensity proportional to norm of density gradi- 
ent). 


2.5 Test computation - flat plate 

A test computation of a turbulent boundary layer along a flat plate at = 3 
has been performed. The Reynolds number is Res 1 = 10000, where 8\ is the lam- 
inar displacement thickness corresponding to the inflow station, which is also the 
reference length. We take as reference length the displacement thickness from a lam- 
inar similarity solution since it is uniquely defined corresponding to a downstream 
station measured from the plate leading edge. The flow parameters are given in 
Table 2. Discretization is N x = 351, N y = 41 and N z = 121. 

The inflow data are generated from the temporal simulation data of Guo and 
Adams (1994) using Taylor’s hypothesis. Initial condition is a laminar similarity 
solution which is also the reference solution used in the sponge region 48 < x < 56 
(Adams, 1994). The computation extends over 8000 time steps. Time step size is 
about A t = 0.1069 t + . The output data are sampled over the final 4400 time steps, 
starting after the inflow plane has been convected through the outflow. The time 
sampling interval is about 470 t+. 

For a comparison we refer to the inflow boundary layer profile of the experimental 
data at higher Reynolds number for a 25° compression corner of Zheltovodov et 
al. (1990). In Table 3 we compare data from simulation and experiment. In 
Fig. 4 we compare mean flow profiles (spanwise and ensemble averaged) at the 
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TABLE 2. Flow parameters for Moo = 3 flat plate. 


quantity 

value 

comment 


115K 


Moo 

3 


Pr 

0.72 


K 

1.4 


R 

287.03 


v’oo 

7.98 • 10 -6 kg/m/s 


s* 

110.4K 


Re»i 

10000 


s? 

4.0830 • 10 -4 m 


£o 

338.32 

inflow station 

L x 

56 

streamwise box-length 

Ly 

4 

spanwise box-length 

L z 

25 

wall-normal box-length 


TABLE 3. Boundary layer data for flat plate, Cf is the skin friction coefficient, 
is the friction velocity, / + is the wall unit, A+ is the grid spacing in wall units 
(for the wall-normal direction z it is the distance of the first point away from the 
wall), and 6i is the turbulent displacement thickness. 


quantity 

x = 10.08 

x = 25.12 

x = 40.00 

exp 

c f 

0.27 • 10" 2 

0.28 • 10“ 2 

0.26 • 10‘ 2 

0.15 • 10~ 2 


0.0595 

0.0605 

0.0576 

0.0442 

7+ 

0.0107 

0.0101 

0.0098 

0.0075 

A+ 

17.41 7+ 

GO 

19.02 1+ 

- 

A+ 

9.33 /+ 

9.89 l + 

10.19 /+ 

- 

A^ 

4.89 1+ 

5.18 /+ 

5.34 7+ 

- 

S i 

1.53 

1.56 

1.55 

1.93 


same streamwise stations as in Table 3 with the experiment. 

In general the quality of the simulation data is unsatisfactory. This is due to 
several reasons. One is the large distance of the first grid point away from the 
wall, which results in a poor approximation of wall-normal gradients. Another 
is the relatively small streamwise extent of the computational domain, which is 
only about 10 turbulent boundary layer thicknesses, considering the fact that the 
outflow sponge affects about another 1.5 boundary layer thicknesses, even less. The 
downstream extent of the inflow transient cannot be clearly assessed. Also, we 
make the same observation as in Guo and Adams (1994) that there is a mass defect 
visible in the profiles from the simulation in the lower boundary layer half. This 
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FIGURE 4. Mean flow profiles for Moo — 3 flat plate. Symbols: , x = 10; 

, x = 25; , x — 40; o , exp. 

is attributed to the too small streamwise extent, which apparently does not allow 
for the appropriate evolution of streamwise streaks. The computational cost was 
32 ns/(N P oint Ntimestep ) for a single CRAY C90 CPU. 

3. Future plans 

From the numerical experiments mentioned in the previous section, we estimate 
a set of parameters where a direct numerical simulation of a compression corner 
is feasible. A direct numerical simulation at these parameters will be attempted 
while an accompanying large-eddy simulation is under consideration by K. Mahesh 

(CTR). 

3. 1 DNS parameters and cost estimate 

The Reynolds number with respect to the turbulent displacement thickness at 
inflow is about 6000. Turbulent boundary layer thickness and turbulent displace- 
ment thickness can be estimated as about 600 / + and 210 / + , respectively. With 
an expected discretization of N x — 601, N y = 51, and N z = 141, we estimate 
A* = 15 / + and A z = 10 l + . With an estimated T pat s = 300 t + for the inflow 
plane to be convected through the domain, a time step of about At = 0.06 t + , and 
a code performance of about Z8ns/(N po i„t Ntimestep ) on a single CRAY C90 CPU, 
we require an estimated 265 hours per T pa> ,. 
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TABLE 4. Flow parameters for = 3 ramp. 


quantity 

value 

comment 

2S, 

115K 


Moo 

3 


Pr 

0.72 


K 

1.4 


R 

287.03 


»*oo 

7.98 • 10 -6 kg/m/s 


S* 

110.4K 


Retx 

4000 



1.6331 • 10 -4 m 

reference length 

& 

135.33 

inflow station 

Lx 

45 

length of first ramp segment 

l 2 

45 

length of second ramp segment 

4 > 

18 

ramp deflection angle 
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Appendix A, Split form of the convective fluxes 

A typical indication of underresolution (thus of aliasing errors) of a direct numer- 
ical simulation of a compressible flow, solving the compressible Navier-Stokes equa- 
tions, is the appearance of regions with negative temperature (or pressure). This 
is related to a local imbalance of internal (potential) and kinetic energy, caused 
mostly by aliasing errors. It has been observed by Blaisdell et aL (1991) that 
for the pseudospectral computation of derivatives of convolutions of dependent 
variables, as d(fg), the aliasing error is reduced by using the identity d(fg ) = 
1/2 (d(fg) + 1/2 + fdg + 1/2 gdf. For finite-difference schemes the coefficients of 
the discrete Fourier series for the derivative have to be multiplied by the integer 
modified wavenumber, which becomes a function of the integer wavenumber; for 
dissipative schemes this modified wavenumber is complex. In this appendix we 
briefly investigate the effect of a split form of the convective fluxes for a dissipative 
finite-difference schemes. From numerical experimentation with coarsely resolved 
computations for a flat plate, we see that for the upwind scheme used above aliasing 
errors are even more critical for the split formulation than for the conservative form. 

First we derive the expressions for pseudospectral convolution in terms of discrete 
Fourier series for a Fourier scheme (in the following the summations £ n +m=* anc ^ 
X)n+m=*±N are always to be taken over m, n = —N/ 2, ..., — N/2 - 1) 


-l 


9 X {f 9) — ^ y J ^ y fm9n + y ^ fm9n 


k~-^- \m+n=k 


m+n=k±N 


) e ikx 


(A.1) 


and 


N y 

E ( E 

k=-^ \m+n=k 


\d x Ug) + \gd*f + \fd*g = 

.1, 


-(k + m + n)f m g„ + ^ i-(k + m + n)f m g„ \ e 

m+n=k±N ’ 


i kx 


£L-\ 


= X] ( E i kf m g n + )/«*$» 


k=~^Y \m+Ti=k 


m+n=k±N 

For a finite difference schemes this reads 


^ e ikx . 


(A.2) 


— -1 

*> A 


dx(fg) = ^ i k(k) f fm9n + Y1 f m 9 n ) ^ (^- 3 ) 

\m+n=k m+n=k±N 


k=-Z- 

K— 2 


') 


\d x (fg) + \d z fg + \fd z g = 


and 
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= J I J + M m ) + fl i T i))fm9 n + 

k= — ^- \m+n=k 

+ *§(*(*) + *™( m ) + fi(n))f m g„ j e' kx : 

m+n=k±N / 


(AA) 


Using the symmetry properties of the dispersion function and of the dissipation 
function, we restrict the following to 0 < £ < 7r. We approximate the modified 
wavenumber £(£) piecewise linearly by 


^ l + , if S < e < tt ’ 

The integer modified wavenumber fe(fc) is then obtained as 

h , x / fc , if fc < A' 

*v fc ) — S N-2fc ft- | : 2A'-2fc T\ \f K < k < — 

l N- 2 K* + 1 N- 2 K U ' 11 A ^ * - 2 

For the split formulation we get: 

Case (1) -K <k< K: 

\d x (fg)+\d x fg+ 1 -f d t g = 

N 
T 


(A.5) 


(A.6) 


- E 

Case (2) A < Jfc < f : 


i ( 2 J^,fmg n + 

E 1 

K)/- 

•'ll 

0n 1 

\m+n=k i 

m+n—k±N 


/. 


JJbx 


(A.7) 


ia,(/s) + + j/a,s = 


= E 


* 2 


E 3 K — 2k p A 

— j\r — 9 A" ^ + E 


3A - — 2k : . a 
- U Jm9n 4" 



+ i 


V- \ N - 2k rsi , , V- * 

) , N — 2 * T * rr ** f™9n 

-v— 

II 


m+n=it > 


„ iV — 2A' 

m+n= fc±N s ■ ■ v ■■■ 


III 


e itr . (A.8) 


For the non-split formulation we get: 
Case (1) —K < k < K: 


d z (fg) = 
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N 

T' 


-1 

= E 

*=-f 

Case (2.) K < k < f : 


E sj^fm9n + ^ kfm<jn 

.m+n=fc iv m+n=fc±N 


a,(/ 5 ) = 




(A.9) 


JV 

T“ 


--1 

- £ 

*=-* 


2A — 2k n ~ A v*^ 2A — 2fc - * 

E " AT IToK D f m9n + E “ N -2K D ^ m9n + 


+ i 


/ \ 

N — 2k Tr ^ -/V — 2 A; - * 

W^2K K mg " + ^ TV - 2K A ^ m9n 

m+n=k s. ^ ✓ m-|-n=HN 


/j 


JJfcx 


(A10) 



•N/2 0 N/2 k 


FIGURE 5. Sketch of the dispersion for non-split and split formulation. 
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Figure 6. Sketch of the dissipation for non-split and split formulation. 
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The factors in the above Fourier sums have character of modified wavenumbers 
and represent the wave properties (dispersion and dissipation) of the particular 
Fourier mode. We now inspect particular factors. Considering the first terms on 
the left-hand sides first, we see that the split formulation generates the spurious 
wave (II, VI) while the non-split formulation generates (V,VIII). From the disper- 
sion shown in Fig. 5 it is evident that II and V contribute by aliasing to the resolved 
spectrum. The spurious waves (II, VI) from the split formulation however, are par- 
tially amplified (negative dissipation) while the spurious waves from the non-split 
formulation are damped. From the second terms on the right-hand side we see that 
the split formulation generates another pair of spurious waves which contribute to 
the resolved spectrum by aliasing which is also amplified (III, VII). We conclude 
that the non-split formulation for an upwind scheme can exhibit spurious waves 
which are amplified contrary to the non-split form. This is due to the fact that the 
modified wavenumber for dissipative schemes is complex. 



